{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thomson Scattering: Spectral Density"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[thomson]: ../../diagnostics/thomson.rst\n",
    "[spectral-density]: ../../api/plasmapy.diagnostics.thomson.spectral_density.rst#spectral-density\n",
    "[sheffield]: https://www.sciencedirect.com/book/9780123748775/plasma-scattering-of-electromagnetic-radiation\n",
    "\n",
    "The [thomson.spectral_density][spectral-density] function calculates the [spectral density function S(k,w)][sheffield], which is one of several terms that determine the scattered power spectrum for the Thomson scattering of a probe laser beam by a plasma. In particular, this function calculates $S(k,w)$ for the case of a plasma consisting of one or more ion species and electron populations under the assumption that all of the ion species and the electron fluid have Maxwellian velocity distribution functions and that the combined plasma is quasi-neutral. In this regime, the spectral density is given by the equation:\n",
    "\n",
    "\\begin{equation}\n",
    "S(k,\\omega) = \\sum_e \\frac{2\\pi}{k} \\bigg |1 - \\frac{\\chi_e}{\\epsilon} \\bigg |^2 f_{e0,e}\\bigg ( \\frac{\\omega}{k} \\bigg ) + \\sum_i \\frac{2\\pi Z_i}{k} \\bigg | \\frac{\\chi_e}{\\epsilon} \\bigg |^2 f_{i0, i} \\bigg ( \\frac{\\omega}{k} \\bigg )\n",
    "\\end{equation}\n",
    "\n",
    "where $\\chi_e$ is the electron component susceptibility of the plasma and $\\epsilon = 1 + \\sum_e \\chi_e + \\sum_i \\chi_i$ is the total plasma dielectric function (with $\\chi_i$ being the ion component of the susceptibility), $Z_i$ is the charge of each ion, $k$ is the scattering wavenumber, $\\omega$ is the scattering frequency, and the functions $f_{e0,e}$ and $f_{i0,i}$ are the Maxwellian velocity distributions for the electrons and ion species respectively.\n",
    "\n",
    "Thomson scattering can be either non-collective (the scattered spectrum is a linear sum of the light scattered by individual particles) or collective (the scattered spectrum is dominated by scattering off of collective plasma waves). The [thomson.spectral_density][spectral-density] function can be used in both cases. These regimes are delineated by the dimensionless constant $\\alpha$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\alpha = \\frac{1}{k \\lambda_{De}}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\lambda_{De}$ is the Debye length. $\\alpha > 1$ corresponds to collective scattering, while $\\alpha < 1$ corresponds to non-collective scattering. Depending on which of these regimes applies, fitting the scattered spectrum can provide the electron (and sometimes ion) density and temperature. Doppler shifting of the spectrum can also provide a measurement of the drift velocity of each plasma species.\n",
    "\n",
    "For a detailed explanation of the underlying physics (and derivations of these expressions), see [\"Plasma Scattering of Electromagnetic Radiation\" by Sheffield et al.][sheffield]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import astropy.units as u\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from plasmapy.diagnostics import thomson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct parameters that define the Thomson diagnostic setup, the probing beam and scattering collection.  These parameters will be used for all examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The probe wavelength can in theory be anything, but in practice integer frequency multiples of the Nd:YAG wavelength\n",
    "# 1064 nm are used (532 corresponds to a frequency-doubled probe beam from such a laser).\n",
    "probe_wavelength = 532 * u.nm\n",
    "\n",
    "# Array of wavelengths over which to calculate the spectral distribution\n",
    "wavelengths = (\n",
    "    np.arange(probe_wavelength.value - 60, probe_wavelength.value + 60, 0.01) * u.nm\n",
    ")\n",
    "\n",
    "# The scattering geometry is defined by unit vectors for the orientation of the probe laser beam (probe_n) and\n",
    "# the path from the scattering volume (where the measurement is made) to the detector (scatter_n).\n",
    "# These can be setup for any experimental geometry.\n",
    "probe_vec = np.array([1, 0, 0])\n",
    "scattering_angle = np.deg2rad(63)\n",
    "scatter_vec = np.array([np.cos(scattering_angle), np.sin(scattering_angle), 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to calculate the scattered spectrum, we must also include some information about the plasma. For this plot we'll allow the ``fract``, ``ion_species``, ``fluid_vel``, and ``ion_vel`` keywords to keep their default values, describing a single-species H+ plasma at rest in the laboratory frame. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx-thumbnail": {
     "tooltip": "Thomson Scattering"
    }
   },
   "outputs": [],
   "source": [
    "ne = 2e17 * u.cm**-3\n",
    "T_e = 12 * u.eV\n",
    "T_i = 10 * u.eV\n",
    "\n",
    "alpha, Skw = thomson.spectral_density(\n",
    "    wavelengths,\n",
    "    probe_wavelength,\n",
    "    ne,\n",
    "    T_e=T_e,\n",
    "    T_i=T_i,\n",
    "    probe_vec=probe_vec,\n",
    "    scatter_vec=scatter_vec,\n",
    ")\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(wavelengths, Skw, lw=2)\n",
    "ax.set_xlim(probe_wavelength.value - 10, probe_wavelength.value + 10)\n",
    "ax.set_ylim(0, 1e-13)\n",
    "ax.set_xlabel(r\"$\\lambda$ (nm)\")\n",
    "ax.set_ylabel(\"S(k,w)\")\n",
    "ax.set_title(\"Thomson Scattering Spectral Density\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example Cases in Different Scattering Regimes\n",
    "\n",
    "We will now consider several example cases in different scattering regimes. In order to facilitate this, we'll set up each example as a dictionary of plasma parameters:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A single-species, stationary hydrogen plasma with a density and temperature that results in a scattering spectrum dominated by scattering off of single electrons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "non_collective = {\n",
    "    \"name\": \"Non-Collective Regime\",\n",
    "    \"n\": 5e15 * u.cm**-3,\n",
    "    \"T_e\": 40 * u.eV,\n",
    "    \"T_i\": np.array([10]) * u.eV,\n",
    "    \"ions\": [\"H+\"],\n",
    "    \"electron_vel\": np.array([[0, 0, 0]]) * u.km / u.s,\n",
    "    \"ion_vel\": np.array([[0, 0, 0]]) * u.km / u.s,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A single-species, stationary hydrogen plasma with a density and temperature that result in weakly collective scattering (scattering parameter $\\alpha$ approaching 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weakly_collective = {\n",
    "    \"name\": \"Weakly Collective Regime\",\n",
    "    \"n\": 2e17 * u.cm**-3,\n",
    "    \"T_e\": 20 * u.eV,\n",
    "    \"T_i\": 10 * u.eV,\n",
    "    \"ions\": [\"H+\"],\n",
    "    \"electron_vel\": np.array([[0, 0, 0]]) * u.km / u.s,\n",
    "    \"ion_vel\": np.array([[0, 0, 0]]) * u.km / u.s,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A single-species, stationary hydrogen plasma with a density and temperature that result in a spectrum dominated by multi-particle scattering, including scattering off of ions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "collective = {\n",
    "    \"name\": \"Collective Regime\",\n",
    "    \"n\": 5e17 * u.cm**-3,\n",
    "    \"T_e\": 10 * u.eV,\n",
    "    \"T_i\": 4 * u.eV,\n",
    "    \"ions\": [\"H+\"],\n",
    "    \"electron_vel\": np.array([[0, 0, 0]]) * u.km / u.s,\n",
    "    \"ion_vel\": np.array([[0, 0, 0]]) * u.km / u.s,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A case identical to the collective example above, except that now the electron fluid has a substantial drift velocity parallel to the probe laser and the ions have a drift (relative to the electrons) at an angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drifts = {\n",
    "    \"name\": \"Drift Velocities\",\n",
    "    \"n\": 5e17 * u.cm**-3,\n",
    "    \"T_e\": 10 * u.eV,\n",
    "    \"T_i\": 10 * u.eV,\n",
    "    \"ions\": [\"H+\"],\n",
    "    \"electron_vel\": np.array([[700, 0, 0]]) * u.km / u.s,\n",
    "    \"ion_vel\": np.array([[-600, -100, 0]]) * u.km / u.s,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A case identical to the collective example, except that now the plasma consists 25% He+1 and 75% C+5, and two electron populations exist with different temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "two_species = {\n",
    "    \"name\": \"Two Ion and Electron Components\",\n",
    "    \"n\": 5e17 * u.cm**-3,\n",
    "    \"T_e\": np.array([50, 10]) * u.eV,\n",
    "    \"T_i\": np.array([10, 50]) * u.eV,\n",
    "    \"efract\": np.array([0.5, 0.5]),\n",
    "    \"ifract\": np.array([0.25, 0.75]),\n",
    "    \"ions\": [\"He-4 1+\", \"C-12 5+\"],\n",
    "    \"electron_vel\": np.array([[0, 0, 0], [0, 0, 0]]) * u.km / u.s,\n",
    "    \"ion_vel\": np.array([[0, 0, 0], [0, 0, 0]]) * u.km / u.s,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "examples = [non_collective, weakly_collective, collective, drifts, two_species]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each example, plot the spectral distribution function over a large range to show the broad electron scattering feature (top row) and a narrow range around the probe wavelength to show the ion scattering feature (bottom row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=len(examples), nrows=2, figsize=[25, 10])\n",
    "fig.subplots_adjust(wspace=0.4, hspace=0.4)\n",
    "\n",
    "lbls = \"abcdefg\"\n",
    "\n",
    "for i, x in enumerate(examples):\n",
    "    alpha, Skw = thomson.spectral_density(\n",
    "        wavelengths,\n",
    "        probe_wavelength,\n",
    "        x[\"n\"],\n",
    "        T_e=x[\"T_e\"],\n",
    "        T_i=x[\"T_i\"],\n",
    "        ifract=x.get(\"ifract\"),\n",
    "        efract=x.get(\"efract\"),\n",
    "        ions=x[\"ions\"],\n",
    "        electron_vel=x[\"electron_vel\"],\n",
    "        ion_vel=x[\"ion_vel\"],\n",
    "        probe_vec=probe_vec,\n",
    "        scatter_vec=scatter_vec,\n",
    "    )\n",
    "\n",
    "    ax[0][i].axvline(x=probe_wavelength.value, color=\"red\")  # Mark the probe wavelength\n",
    "    ax[0][i].plot(wavelengths, Skw)\n",
    "    ax[0][i].set_xlim(probe_wavelength.value - 15, probe_wavelength.value + 15)\n",
    "    ax[0][i].set_ylim(0, 1e-13)\n",
    "    ax[0][i].set_xlabel(r\"$\\lambda$ (nm)\")\n",
    "\n",
    "    ax[0][i].set_title(lbls[i] + \") \" + x[\"name\"] + f\"\\n$\\\\alpha$={alpha:.4f}\")\n",
    "\n",
    "    ax[1][i].axvline(x=probe_wavelength.value, color=\"red\")  # Mark the probe wavelength\n",
    "    ax[1][i].plot(wavelengths, Skw)\n",
    "    ax[1][i].set_xlim(probe_wavelength.value - 1, probe_wavelength.value + 1)\n",
    "    ax[1][i].set_ylim(0, 1.1 * np.max(Skw.value))\n",
    "    ax[1][i].set_xlabel(r\"$\\lambda$ (nm)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plots of the spectral density function (Skw) which determines the amount of light scattered into different wavelengths.\n",
    "\n",
    "a. In the non-collective regime only the electron feature is visible.\n",
    "\n",
    "b. In the weakly collective regime (alpha approaches 1) an ion feature starts to appear and the electron feature is distorted\n",
    "\n",
    "c. In the collective regime both features split into two peaks, corresponding to scattering off of forward and backwards propagating plasma oscillations. \n",
    "\n",
    "d. The introduction of drift velocities introduces several Doppler shifts in the calculations, resulting in a shifted spectrum.\n",
    "\n",
    "e. Including multiple ion and electron populations modifies the ion and electron features respectively.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
